home *** CD-ROM | disk | FTP | other *** search
/ HAKERIS 11 / HAKERIS 11.ISO / linux / system / LinuxConsole 0.4 / linuxconsole0.4install-en.iso / guile0.4.lcm / share / guile / slib / minimize.scm < prev    next >
Encoding:
Text File  |  2004-01-06  |  3.9 KB  |  115 lines
  1. ;;; "minimize.scm" finds minimum f(x) for x0 <= x <= x1.
  2. ;;; Author: Lars Arvestad
  3. ;;;
  4. ;;; This code is in the public domain.
  5.  
  6. ;;@noindent
  7. ;;
  8. ;;The Golden Section Search
  9. ;;@footnote{David Kahaner, Cleve Moler, and Stephen Nash
  10. ;;@cite{Numerical Methods and Software}
  11. ;;Prentice-Hall, 1989, ISBN 0-13-627258-4}
  12. ;;algorithm finds minima of functions which
  13. ;;are expensive to compute or for which derivatives are not available.
  14. ;;Although optimum for the general case, convergence is slow,
  15. ;;requiring nearly 100 iterations for the example (x^3-2x-5).
  16. ;;
  17. ;;@noindent
  18. ;;
  19. ;;If the derivative is available, Newton-Raphson is probably a better
  20. ;;choice.  If the function is inexpensive to compute, consider
  21. ;;approximating the derivative.
  22.  
  23. ;;@body
  24. ;;
  25. ;;@var{x_0} are @var{x_1} real numbers.  The (single argument)
  26. ;;procedure @var{f} is unimodal over the open interval (@var{x_0},
  27. ;;@var{x_1}).  That is, there is exactly one point in the interval for
  28. ;;which the derivative of @var{f} is zero.
  29. ;;
  30. ;;@0 returns a pair (@var{x} . @var{f}(@var{x})) where @var{f}(@var{x})
  31. ;;is the minimum.  The @var{prec} parameter is the stop criterion.  If
  32. ;;@var{prec} is a positive number, then the iteration continues until
  33. ;;@var{x} is within @var{prec} from the true value.  If @var{prec} is
  34. ;;a negative integer, then the procedure will iterate @var{-prec}
  35. ;;times or until convergence.  If @var{prec} is a procedure of seven
  36. ;;arguments, @var{x0}, @var{x1}, @var{a}, @var{b}, @var{fa}, @var{fb},
  37. ;;and @var{count}, then the iterations will stop when the procedure
  38. ;;returns @code{#t}.
  39. ;;
  40. ;;Analytically, the minimum of x^3-2x-5 is 0.816497.
  41. ;;@example
  42. ;;(define func (lambda (x) (+ (* x (+ (* x x) -2)) -5)))
  43. ;;(golden-section-search func 0 1 (/ 10000))
  44. ;;      ==> (816.4883855245578e-3 . -6.0886621077391165)
  45. ;;(golden-section-search func 0 1 -5)
  46. ;;      ==> (819.6601125010515e-3 . -6.088637561916407)
  47. ;;(golden-section-search func 0 1
  48. ;;                       (lambda (a b c d e f g ) (= g 500)))
  49. ;;      ==> (816.4965933140557e-3 . -6.088662107903635)
  50. ;;@end example
  51.  
  52. (define golden-section-search
  53.   (let ((gss 'golden-section-search:)
  54.     (r (/ (- (sqrt 5) 1) 2)))    ; 1 / golden-section
  55.     (lambda (f x0 x1 prec)
  56.       (cond ((not (procedure? f)) (slib:error gss 'procedure? f))
  57.         ((not (number? x0)) (slib:error gss 'number? x0))
  58.         ((not (number? x1)) (slib:error gss 'number? x1))
  59.         ((>= x0 x1) (slib:error gss x0 'not '< x1)))
  60.       (let ((stop?
  61.          (cond
  62.           ((procedure? prec) prec)
  63.           ((number? prec)
  64.            (if (>= prec 0)
  65.            (lambda (x0 x1 a b fa fb count) (<= (abs (- x1 x0)) prec))
  66.            (if (integer? prec)
  67.                (lambda (x0 x1 a b fa fb count) (>= count (- prec)))
  68.                (slib:error gss 'integer? prec))))
  69.           (else (slib:error gss 'procedure? prec))))
  70.         (a0 (+ x0 (* (- x1 x0) (- 1 r))))
  71.         (b0 (+ x0 (* (- x1 x0) r)))
  72.         (delta #f)
  73.         (fmax #f)
  74.         (fmin #f))
  75.     (let loop ((left x0)
  76.            (right x1)
  77.            (a a0)
  78.            (b b0)
  79.            (fa (f a0))
  80.            (fb (f b0))
  81.            (count 1))
  82.       (define finish
  83.         (lambda (x fx)
  84.           (if (> fx fmin) (slib:warn gss fx 'not 'min (list '> fmin)))
  85.           (if (and (> count 9) (or (eqv? x0 left) (eqv? x1 right)))
  86.           (slib:warn gss 'min 'not 'found))
  87.           (cons x fx)))
  88.       (case count
  89.         ((1)
  90.          (set! fmax (max fa fb))
  91.          (set! fmin (min fa fb)))
  92.         ((2)
  93.          (set! fmin (min fmin fa fb))
  94.          (if (eqv? fmax fa fb) (slib:error gss 'flat? fmax)))
  95.         (else
  96.          (set! fmin (min fmin fa fb))))
  97.       (cond ((stop? left right a b fa fb count)
  98.          (if (< fa fb)
  99.              (finish a fa)
  100.              (finish b fb)))
  101.         ((< fa fb)
  102.          (let ((a-next (+ left (* (- b left) (- 1 r)))))
  103.            (cond ((and delta (< delta (- b a)))
  104.               (finish a fa))
  105.              (else (set! delta (- b a))
  106.                    (loop left b a-next a (f a-next) fa
  107.                      (+ 1 count))))))
  108.         (else
  109.          (let ((b-next (+ a (* (- right a) r))))
  110.            (cond ((and delta (< delta (- b a)))
  111.               (finish b fb))
  112.              (else (set! delta (- b a))
  113.                    (loop a right b b-next fb (f b-next)
  114.                      (+ 1 count))))))))))))
  115.